Infinite size density matrix renormalization group, revisited 
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I revisit the infinite-size variant of the Density Matrix Renormalization Group (iDMRG) algorithm 
for obtaining a fixed-point translationally invariant matrix product wavefunction in the context of 
one-dimensional quantum systems. A crucial ingredient of this algorithm is an efficient transfor- 
mation for obtairung the matrix elements of the wavefunction as the lattice size is increased, and 
I introduce here a versatile transformation that is demonstrated to be much more effective than 
previous versions. The resulting algorithm has a surprisingly close relationship to Vidal's Time 
Evolving Block Decimation for infinite systems, but allows much faster convergence. Access to a 
translationally invariant matrix product state allows the calculation of correlation functions based 
on the transfer matrix, which directly gives the spectrum of all correlation lengths. I also show some 
advantages of the Matrix Product Operator (MPO) technique for constructing expectation values 
of higher moments, such as the exact variance {{H — E)'^). 
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I. INTRODUCTION 

For many years, the density-matrix renormalization- 
group algorithniiiS, (DMRG) has been used in an in- 
creasing diversity of applications in fields such as con- 
densed matter physics, condensed matter physics, quan- 
tum chemistry, nuclear physics and quantum information 
science. The DMRG algorithm is used to obtain a varia- 
tional approximation for groundstates and low-lying ex- 
cited states, with later variants being developed for finite- 
temperature and dynamical properties. The basic algo- 
rithm comes in two varieties. The most commonly used 
algorithm is the so-called finite-size algorithm, whereby 
a matrix-product state (MPS) wavefunction on a finite- 
size lattice (defined through a distinct set of matrices 
A*", for each lattice site n — 1,2,...,L, where |s„) 
is a d-dimensional basis for the local Hilbert space at 
site n) is iteratively improved until convergence. The 
infinite-size algorithm, on the other hand, grows the lat- 
tice by one or more sites each iteration and produces, 
at the fixed point, a wavefunction that is translation- 
ally invariant (possibly with a non-trivial unit cell) . This 
approach is useful for probing quantities in the thermo- 
dynamic limit, without the influence of boundary effects, 
and such a translationally invariant state can readily be 
used as an initial state to study real-time evolution of a 
homogeneous quench as well as local (translation symme- 
try breaking) perturbations. Recently, interest in MPS 
algorithms that obtain a thermodynamic fixed point has 
been rekindled with the invention of a variant of Vidal's 
Time-Evolving Block Decimation (TEBD) algorithniii 
that applies directly to an infinite size wavefunction to 
give iTEBD^^, the resulting wavefunction being trans- 
lationally invariant by construction (although with the 
restriction that the unit cell must be at least two sites in 
size). The iTEBD utilizes imaginary time evolution via a 
Trotter-Suzuki decomposition^, amounting to using the 
power method^^ to obtain the groundstate eigenvector of 
an approximation of the exponential of the Hamiltonian. 



However the power method is not an efficient eigensolver, 
and the use of the Trotter-Suzuki decomposition requires 
careful scaling of the time-step to zero in order to obtain 
a properly converged state. Hence, all else being equal, 
a DMRG approach where an efficient local eigensolver 
is used to find a variationally optimal state ought to be 
much more robust and efficient. 

In the past, the main use for the infinite size DMRG 
algorithm however has not been to achieve a fixed point, 
but rather to grow an initial wavefunction of a particular 
size as an initialization step for the finite-size algorithm^i. 
As an algorithm for obtaining a translationally invariant 
fixed point, infinite-size DMRG (hereafter referred to as 
iDMRG) has been virtually abandoned. There are sev- 
eral reasons for this: Firstly, finite-size DMRG (fDMRG) 
readily allows the application finite-size scaling"^^ with 
respect to the lattice size, whereas for iDMRG the only 
available parameter is the number of states kept m. In 
fact, this is not a problem because the scaling relations in 
terms of the number of states have already been investi- 
gated by several authors2L22ii^, although having received 
insufficient attention this approach is not yet widely uti- 
lized. In fact, having only a single scaling parameter 
is a marked advantage of iDMRG over finite size ap- 
proaches where either the number of states kept must 
be sufficiently large that the truncation effects are negli- 
gible, which may be impractical, or m must be considered 
as an additional parameter and scaling performed with 
respect to both parameters, which in the past has not 
always been done correctly (see Ref. 27 for the correct 
approach) . 

A second disadvantage of iDMRG over fDMRG is a 
technical limitation of the algorithm, namely the lack 
of a good initial guess vector for the local eigensolver 
at each DMRG iteration. In contrast, the development 
of the wavefunction acceleration procedure in fDMRG 
to transform the basis for the superblock wavefunction 
from one iteration to be the initial guess vector for the 
next iteration, was a breakthrough in the development of 
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finite-size DMRG that improves the efficiency of the al- 
gorithm by orders of magnitude. This transformation is 
essentially implicit in the matrix product formulation of 
DMRG^iSiil. From this point of view, DMRG is an algo- 
rithm for locally updating the matrix elements of an MPS 
one (or a few) lattice sites at a time, and the wavefunc- 
tion transformation results from the orthonormalization 
the basis which is necessary to ensure numerical stability. 
However, this transformation only applies to finite-size 
MPS, and the apparent lack of a similar transformation 
for iDMRG has been the focus of several studies in the 
last decade or more. 

The problem that needs to be solved to obtain a wave- 
function transformation for iDMRG is to find a proce- 
dure for predicting the sign of the new basis elements 
when the lattice size is increased. This sign structure 
has two components, being the intrinsic nodal-structure 
of the wavefunction, as well as a phase ambiguity of the 
eigenvectors of the density matrix (or singular value de- 
composition). There have been several works that at- 
tempt to obtain a transformation algorithm, either us- 
ing Marshall's sign rule^^ to predict the sign structure of 
the wavefunctiorti^, or directly comparing the sign struc- 
ture from different step a^^i^° . These approaches have had 
mixed success, in that they rely on a 1 — 1 matching of 
density matrix eigenvalues (equivalently, singular values 
arising from the Schmidt coefficients) from successive it- 
erations and attempting to match signs of the eigenvec- 
tors (say, by choosing that the first component of the vec- 
tor will always be positive). Where applicable, the Mar- 
shall's sign rule approach is very effectivei^, but the need 
to know the nodal structure of the wavefunction a priori 
is a severe limitation, eg. this approach breaks down com- 
pletely for frustrated systems. The approach of Qin and 
Lou does lead to a wavefunction that is somewhat bet- 
ter than a random state, and in this sense this approach 
is useful in reducing the number of iterations required 
by the local eigensolver. However, repairing eigenvector 
signs in this manner is not robust and this procedure 
has not been demonstrated to be effective in obtaining 
a fixed point wavefunction; that is, in the studies utiliz- 
ing this approacbi2i2£ the algorithm does not converge 
to a point where the transformation reproduces the fixed 
point wavefunction with full fidelity. 

It may be a surprise, therefore, that a solution to most 
of these problems exists in the literature and has done for 
quite some time, pre-dating even the wavefunction trans- 
formation in fDMRG, under the name of Product Wave- 
Function Renormalization Group (PWFRG)^i2ii25j26^ 
In this approach, a recurrence relation is set up to obtain 
the translation operator that shifts the wavefunction by 
one site (or unit cell, which may be an arbitrary number 
of lattice sites) . Using this operator the initial wavefunc- 
tion for the next iteration can be obtained. However, 
while this procedure is effective in obtaining a transla- 
tionally invariant fixed point, it still relies on matching 
singular values from successive iterations and the rate of 
convergence away from the fixed point is quite poor. This 



is because the wavefunction transformation does not have 
full fidelity until the translation operator itself has con- 
verged, however the translation operator obtained from 
the recurrence relation converges is quite slowly. A main 
result of this paper is the introduction of a new transfor- 
mation that reaches the same fixed point as PWFRG, but 
converges much faster when away from the fixed point, ie. 
when the lattice size is small. In fact, the convergence is 
effective enough that this approach is useful even for ex- 
act diagonalization of small systems. Also in this paper, 
I investigate some approaches to calculating correlation 
functions using the transfer operator. This approach has 
been known for quite some time^i^i^, but has not been 
widely utilized. Thirdly, the more recent development of 
the Matrix Product Operator (MPO) formulation^^, for 
obtaining the exact matrix product representation of the 
Hamiltonian and observables is demonstrated to have sig- 
nificant advantages and I present algorithms for obtain- 
ing expectation values in the infinite-size limit directly 
from the MPO representation. 

The layout of this paper is as follows: In section [Hi 
I describe the basic matrix product approach and no- 
tation, for both wavefunctions and operators. In sec- 
tion mil I describe the iDMRG algorithm itself and some 
discussion of the convergence properties. Section IIVI 
treats the calculation of observables for many cases in- 
cluding local quantities, correlation functions including 
the spectrum of correlation lengths, groundstate fideli- 
ties, and translationally-invariant operators such as the 
energy and the variance. Finally, section |V] contains a 
summary and some concluding remarks. 



II. MATRIX PRODUCT STATES 

We denote an MPS on an i-site lattice by the form, 
Tr ^ ■ • ■ \si)®\s2)(i>---®\sL) . (1) 

The local index Si represents an element of the d- 
dimensional local Hilbert space at site z, and the A^' are 
TO X TO matrices, to being the dimension of the matrix ba- 
sis (this quantity is also often called D, and sometimes 
x). As a short-hand notation, it is convenient to omit 
the basis states | Si) and the summation from Eq. ([T]) and 
speak of an MPS as being defined purely as a sequence 
of matrices. This notation is quite consistent if one re- 
gards the local index Sj as being a ket vector. Thus, 
in shorthand notation, Eq. ([T]) is equivalent to simply 
A'^A'^ • • • A"^. 

In practice, a MPS state with no particular constraints 
on the form of the A-matrices is numerically difficult to 
handle and it is usual to impose some orthonormalization 
constraints (see eg. Ref. [ij for details) , 

^A'^U'^^1, (2) 
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for the left-orthonormalization constraint, and conversely 



J2 B'B'^I ^ 1 



(3) 



for the right-orthonormahty constraint. In this paper, I 
use always to denote a set of matrices satisfying the 
left-orthonormalization constraint Eq. ^ , and to de- 
note a set of matrices satisfying the right constraint Eq. 
([3]). The infinite-size DMRG algorithm is based around 
adding sites to the center of the chain. To this end, it is 
useful to orthonormalize the sites of the MPS from the 
edges, working in towards the center, to give a wavefunc- 
tion of the form, 



A' 



.B' 



(4) 



Here A is a matrix that is left over after all of the sites 
in the MPS have been orthonormalized, which I refer to 
as the center matrix, and this matrix is also equivalent 
to what is called the superblock wavefunction in DMRG 
terminology.^ In Ref. JJi,, I used C for this matrix, here I 
use A to stand for the matrix of singular values (Ai) as 
start towards harmonizing the various notations associ- 
ated with DMRG and MPS approaches. 

These orthonormalization conditions correspond to a 
basis in which the reduced density matrix for the left 
(right) half of the system is diagonal, and are imposed for 
example by performing a singular value decomposition 
treating A^ as a single dm x m matrix (or as a m x dm 
matrix) , 



A^ 
A' 

B = [B^B^ 



(5) 



.B" 



and then performing the singular value decomposition 
A = UAV'^, where A is a m x m diagonal matrix 
such that A^ is the reduced density matrix, and U, V 
are row-unitary matrices. Upon rewriting U back as a 
set of m X m matrices analogously to Eq. ([5|), one ob- 
tains C/" as satisfying the orthogonality constraint Eq. 
(|2|, with a remainder term AV'' appearing on the right 
(which we can, for example, incorporate into the neigh- 
boring A-matrix) . The procedure for obtaining the right- 
orthonormalization constraint is similar, now ending up 
with a remainder matrix appearing on the left. In this 
fashion, the boundary of the left/right orthonormalized 
states can be moved to an arbitrary position in the lat- 
tice, a procedure I refer to as a rotation of the center 
matrix to the left or right. 

The infinite-size DMRG produces is a wavefunction 
that grows by adding sites in the center, which, for the 
example case of a 2-site unit cell, gives a wavefunction at 
the n*™ iteration that is 2n sites long, 



Sites are added until the the A'' and B" matrices in the 
center of the lattice converge to a fixed point. At this 
fixed point, the MPS will be translationally invariant, 
possibly with a non-trivial unit cell, such that the A'^, B^ 
matrices are independent of position. In TEBD the as- 
sociated MPS is always of this translationally invariant 
form, but in iDMRG we will see that the translational 
invariance arises from the wavefunction transformation 
converging to the identity operator at the fixed point. 

The MPS form Eq. ([T|) is related to the canonical form 
used by Vidalii, 



(7) 



by A^ = AF* and B" = F^A*. For the canonical form, 
the orthogonality constraints take the form 



(8) 



where — AA^ and pn — A^A are the reduced den- 
sity matrices for the left and right semi-infinite strips 
respectively. This shows that the F'* and A matrices for 
a canonical MPS are not independent. In terms of the 
A'',B'' matrices, the corresponding fixed point relations 
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(9) 



(6) 



States that have a non-zero quantum number can be 
readily obtained in this scheme, as long as a repeated 
unit cell can be properly defined. For example, for fi- 
nite doping or finite magnetization with a density p/q 
for coprime integers p, q, we use a unit cell of q sites and 
insert a shift in the quantum numbers of —p per unit cell. 
The effect of this shift in the wavefunction (equivalent to 
an auxiliary particle of quantum number — p) is to force 
the physical system to have an average quantum num- 
ber of +p throughout the unit cell, while the entire MPS 
(physical system plus the quantum number shift) has an 
average quantum number of zero and thereby allowing 
a repeated unit cell. This shift can be implemented in 
a variety of equivalent ways, for example by introducing 
explicitly an additional matrix Q into the MPS unit cell 
that has matrix elements {j' \ Q\j) = Sj'j-p, or by shift- 
ing the quantum numbers in the local basis of one site in 
the unit cell by —p. Note that this procedure does not 
enforce that there will be exactly p particles per unit cell, 
which would be an undesirable restriction, but arbitrarily 
large fluctuations are allowed within the constraint that 
the average quantum number on the infinite system is 

p/q- 



A. Matrix product operators 

Ref. ^3 introduces a representation for the Hamilto- 
nian and other associated operators as an exact matrix 
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product operator (MPO). These operators come in sev- 
eral varieties, the main form that we consider here is op- 
erators that are characterized by an MPO representation 
that can be put in a lower (upper) triangular form, and 
therefore represent operators that are polynomial com- 
binations of local operators, hence the expectation value 
of such an operator is a polynomial function of the lat- 
tice size, and in typical cases will be linear, for exam- 
ple giving the energy per site. A second form of MPO 
arises from products of local operators, such as occurs 
in the exponential of a triangular operator. This form 
has an expectation value that is an exponential function 
of the lattice size, hence I call this class exponential 
MPO's. Crosswhite and Bacon^i have also considered 
MPO's from the point of view of a correspondence with 
weighted finite-state automata. 

Firstly, we briefly examine some properties of triangu- 
lar MPO's. The notation we use here is that an MPO is 
a matrix of local operators of dimension M x M. This 
form, after becoming familiar with the notation, is a very 
convenient way of representing an MPO, see Ref. [l3| for 
full details. For simplicity we take here the unit cell to 
be a single site, so that the infinite MPO is represented 
by the same matrix repeated on every site. For cases 
where the unit cell of the Hamiltonian is non-trivial, for 
example an iV-leg ladder, the associated MPO has the 
same periodicity which leads to the technicality of treat- 
ing MPO matrices that may be non-square (trapezoidal, 
in the case of triangular MPO's), but this presents no 
serious difficulty and all of the procedures described here 
carry over to this case straightforwardly. Starting from 
the most general lower-triangular 2x2 MPO, 

^={b'a)^ (10) 

this represents the infinite sum of all terms of the form 

A® ■ • ■ (g) A (g) B (g) C (g) ■ • ■ «) C . (11) 

Note that Ref. 13 used M to denote the MPO matrix, 
whereas here I use W to avoid confusion with M, the 
dimension of the W matrix. By appropriate choices of 
local operators A, B, C, this form can be used for many 
variants of one-particle operators, including bosonic and 
fermionic particle operators at arbitrary lattice momenta, 
or the sum of local interactions. Moving to 3 x 3 matrices, 
we can generalize this form to a bond operator, or by 
including a term on the diagonal, 

(E 0\ 
W=\d C , (12) 
\0 B A) 

which represents the sum of all string-like terms of the 
form 

A®---®A®B®C®---®C®D®E®---®E. (13) 

Thus long-range string operators can be represented just 
as easily as short-range terms. Taking the string opera- 
tor C to be proportional to the identity operator shows 



that non-local but exponentially decaying interactions 
can be used in DMRG with negligible loss of efficiency 
over purely local interactions, a fact that is probably well 
known to DMRG practitioners although I am not aware 
of any studies that make use of this. 

For a concrete example of a triangular MPO, combin- 
ing a bond term and a local term a simple but non-trivial 
Hamiltonian operator is the Ising model in a transverse 
field (ITF), 

^= E '^^'^I + ^E^' (14) 

<i,j> i 

which has the lower-triangular MPO form as 3 x 3 ma- 
trices, 

/ 1 0\ 

\ a'- 0. (15) 
\Xa'= (j^ l) 

The use of MPO's gives a convenient way to construct 
a software code that has no explicit dependence on the 
particular Hamiltonian. The main advantage, however, is 
that the MPO form allows arithmetic operations, in par- 
ticular sums and products of operators are represented 
simply by taking the matrix direct sum and direct prod- 
uct respectively of the MPO representation, as described 
in Ref. [l^. If both operators have a lower-triangular 
MPO representation, then the sum or product is also 
a lower-triangular MPO. This means that it is very easy 
to obtain operators that would take considerable effort to 
construct using the ad-hoc methods traditionally used in 
DMRG and MPS approaches. One such example is the 
variance, {H — E)"^, which replaces the discarded weight 
(the truncation error, in DMRG terminology) as the con- 
vergence measure of choice as it gives a measure that is 
both independent of the details of the algorithm and a 
more reliable indicator of convergence. 

Not all useful operators have the form of a triangu- 
lar MPO. The triangular MPO's represent sums of ei- 
ther short-range or string terms, whereas other kinds 
of operators of interest comprise instead products of 
terms. An example of this is the MPO representation 
of the Trotter-Suzuki approximation of an exponential, 
for example the first order decomposition of the expo- 
nential of a Hamiltonian operator with only nearest- 
neighbor terms into terms acting on even and odd bonds, 
exp[-(3H] ~ exp[-/37?evon] exp[-/3ffodd]- The MPO rep- 
resentation of exp[— /3i?even] has a 2-site unit cell with 
alternating 1 x M and M x 1-dimensional M^-matrices 
representing the product of all of the exponentials for 
the even bonds, and similarly exp[— /?_ffodd] alternates 
M X 1 and 1 x Af matrices. The product of the two 
operators therefore is an M x M MPO that will not be 
lower-triangular. This operator falls into the class of ex- 
ponential MPO's, with an expectation value of the form 
for L sites, where A is the eigenvalue corresponding 
to the fixed point i?" eigenmatrices (see below). 

The MPO formulation gives a natural form for obtain- 
ing the actual matrix elements of an operator in the form 
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required to calculate an observable using an MPS wave- 
function. This is done by obtaining matrices E" and 
that represent the operator matrix elements on the 
left and right halves of the system respectively. Later we 
will see how to obtain these matrices for infinite states di- 
rectly, however let us first consider how to construct these 
matrices for a finite, but arbitrary size system, which pre- 
cisely what is needed for the iDMRG algorithm. Let us 
take the wavefunction of size 2n sites, given by Eq. ([6]). 
The initial , F'^ matrices for the boundaries are taken 
to be 

^0 ~ O/Sa-'right ■ 

where Iidt and /right are the identity operators acting on 
the left- and right-hand edges of the wavefunction. For 
an open boundary condition state, which is always what 
we use in iDMRG, these operators will simply be 1 x 1 
identity matrices. The E°',F^ for the remainder of the 
chain are then defined recursively, 

E^' = Es',s,c.W::,iA^'-^E-_,A^'^ 

Once we have obtained the matrices E" , F^ for the cen- 
ter of the chain, we can perform many operations. For 
example, the action of the operator on the center matrix 
in its projected Hilbert space is 

M^(A„)-^KA„F„^t, (18) 

a 

and the expectation value is given by the matrix inner 
product, 

( A„ I I A„ ) = Tr K A„F,f tAt . (19) 

Note that the action on the center matrix Eq. (fTSl) 
occurs in the projected Hilbert space of the m x m- 
dimensional basis that is the support of A„, and Eq. ITSl) 
in no way represents the action of the operator on the 
entire state (but that action is easy to obtain, see for ex- 
ample Ref . [l3l) . It would be a serious error to expect, for 
example, that the action of the product of two operators 
W1W2 on a wavefunction can be obtained by applying 
separately the E, F matrices for W2 and then the E, F 
matrices for Wi. That action would correspond instead 
to the operator W1PW2, where P is the projector onto 
the m X TO-dimensional Hilbert space. Nevertheless, for 
a properly constructed set of E, F matrices, the expec- 
tation value Eq. (jl9p is exact, because it results from 
the inner product of the action of the operator with the 
original state, for which the fact that this operation is 
being performed within a projected Hilbert space does 
not matter. 

Having seen how to obtain matrix elements in the pro- 
jected Hilbert space of a A matrix in the MPS, the final 
ingredient before presenting the iDMRG algorithm is to 



see how the MPS itself is iteratively updated. For this 
task, we need the action of an operator on not just a 
A matrix, but on one (or more) of the A'', B" them- 
selves. This is most easily demonstrated for the action 
of an MPO on a single site and incorporating the center 
matrix. Thus, taking the wavefunction of Eq. ^ and 

leaving all matrices fixed except and A„, we obtain 
the action of an MPO on these matrices, 

where we have used the explicit matrix elements of the 
W matrix defining the MPO W^^" , with cr, s'^ indices 
running over the d-dimensional local basis and a, /3 in- 
dices running over the M-dimensional matrix basis of 
the MPO. For modifying two or more sites at a time, 
analogous formulas apply with one W matrix appearing 
per site. 

III. IDMRG 

In this section, I describe the iDMRG algorithm and 
discuss some convergence and stability issues. The al- 
gorithm itself is presented in figure [TJ For simplicity I 
have used a unit cell of two sites, however the unit cell 
can be any size, including a single site. If the system is 
parity-symmetric then step 3 can be avoided with the B'^ 
matrices taken to be the transpose of A^ (if good quan- 
tum numbers are used, the quantum numbers in the ma- 
trix basis must be flipped as well, otherwise one ends up 
with a conjugate tensor rather than a normal tensor^'^), 
being careful to use an eigenvalue decomposition rather 
than a singular value decomposition to ensure symmetry 
of the left and right basis vectors. For unit cells larger 
than two sites, the only modification is steps 2 and 3 ro- 
tate the center matrix through the remainder of the unit 
cell, optimizing the matrices comprising the unit cell in 
the process. Thus, this transformation gives a recurrence 
relation whereby a wavefunction on an rt-site lattice and 
a wavefunction on a larger n -I- fc-site lattice for some ar- 
bitrary A: > 0, are used to obtain an approximate wave- 
function for an n -I- 2A:-site lattice. Note that steps 2 and 
3 act independently on and could be performed in 
parallel. 

The key component of this algorithm is step 4, where 
the wavefunction center matrix A*'^;^ is obtained from 
those of the previous two iterations, 

A^tl = A,?A,;^A,^ . (21) 

Note that A^ and A^ are center matrices of the same 
wavefunction |*n), but located at different partitions, 
at the left and right boundaries of the unit cell respec- 
tively. Eq. (j2ip is a straightforward computation that, 
unlike the PWFRG transformation^'^ and sign-rule based 
approaches, can be applied directly without any kind of 
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(initialization) Obtain the wavefunction for a 

two-site lattice | "to) ~ -A'q'AoBq" and a four-site 

lattice 1*1) = Ag^A'l^AiBt'Bo", and set n = 1 
for the start of iterations. 

Rotate the center matrix of | '^n) one step to the 

left, to obtain . . . A^B^"^j^ B^" .... 

Rotate the center matrix of | ^l/n) one step to the 

right, to obtain . . . A^"^^']^+^A^ .... 

(new wavefunction) The trial wavefunction for a 

lattice with size increased by two sites is given by 



I ^n + l) — • • ■ ^n+l ^^n^^n-l^^n-On 



1+1 

n 



5. (eigensolver) Use | "i/n+i) as the initial guess vec- 
tor for an eigensolver, updating the matrices 

v4^+VA„+iB^+"\' to obtain the final 

6. (truncation) Truncate the dimension of A„+i as 
desired, using eg. a singular value decomposition 
or eigenvalue decomposition, and normalize. 

7. (stopping criteria) If A„+i is sufficiently close to 
A^'^, then the fixed point has been reached. Oth- 
erwise increase n ^ n + 1 and return to step 2. 

FIG. 1: The iDMRG algorithm, including the wavefunction 
transformation. 



matching of singular values from different basis sets. This 
is because the basis for the right hand side of is the 
same as the right hand side of A„_i, and similarly the 
left hand side of matches the left hand side of A„_i. 
Thus there is not even any requirement here to diagonal- 
ize the A matrices, they need not even be square. For 
simplicity of calculating the inverse A~l]^ however, one 
can choose a basis in which A„_i is diagonal, which is 
normally obtained as part of the truncation step anyway. 

1 - loverlapl 

s=l/2 heisenberg model, 20 states kept, 2-site inoving average 




200 300 
iteration 

FIG. 2: Fidelity 1 - (*Ji"'''|>I'„) of the wavefunction transfor- 
mation for iDMRG and PWFRG. For comparison, the trun- 
cation error is also included, which sets a bound for the maxi- 
mum possible fidelity. To remove small period two oscillations 
caused by the two-site unit cell, a two-site moving average was 
used for all quantities. 



Figure [2] shows the effectiveness of the iDMRG wave- 
function transformation compared with PWFRG for a 
simple example. This figure shows the fidelity of the 
initial wavefunction with respect to the variational min- 
imum, 1 — for the spin-1/2 Heisenberg model 
with a 2-site unit cell and rn ~ 20 states kept in an 
S'C/(2)-symmetric basi o^^i^^ . The correlation length for 
this basis size is ^ 48.5 lattice sites. The exact ground- 
state of this Hamiltonian is critical, however for a finite 
basis size the correlation length of a matrix product wave- 
function will always be finitCi^iS, The signature of critical 
behavior is a scaling law for the correlation length as the 
basis size is increased, of the form ^ = m", where the 
exponent k ~ 1.3 for the Heisenberg mode L^^'^^ Also in- 
cluded in the figure is a variant of Eq. (PTjl that avoids 
the matrix inverse by taking instead the square root of 
the singular values of A^A^ in the basis where A„_i is 
diagonal, but this transformation is much less effective 
for small systems. As expected, the PWFRG transfor- 
mation converges only for system sizes much greater than 
the correlation length due to the slow convergence of the 
translation operator. 

An interesting feature of Fig. [2] is the excellent perfor- 
mance of the iDMRG transformation for systems much 
smaller than the correlation length where the transfor- 
mation captures the effects of the open boundary condi- 
tions extremely well, even for systems small enough to 
be exactly diagonalizable. But for an exact diagonaliza- 
tion, periodic boundary conditions are more natural. So, 
in Fig. [31 I compare the effectiveness of the transforma- 
tion also for periodic boundary conditions, increasing the 
number of states kept to the limit of the numerical ac- 
curacy of the wavefunction transformation. Using just 
10 Lanczos iterations at each lattice size, we obtain an 
energy for a 30-site system with periodic boundary condi- 
tions of -E = —13.321963058(9), accurate to 10 significant 
figures, showing that this transformation could be useful 
for accelerating exact diagonalization studies. 
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FIG. 3: Fidelity 1 - (\['J,"'''|\E'„) of the wavefunction transfor- 
mation for open and periodic boundary conditions. 

Clearly this algorithm has a similar function to 
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iTEBD in the sense that both algorithms result in 
a translationally-invariant variational approximation to 
the groundstate wavefunction. But in fact the relation- 
ship between iDMRG and iTEBD is much deeper. If we 
replace step 5 by an alternative, 

5. (evolution) Obtain | 'I'n+i) by applying a sin- 
gle bond evolution operator to sites s^+i and (22) 

sn+i of 

then the iTEBD algorithm is recovered exactly, where the 
translationally invariant wavefunction at the n**^ step is 
given, in notation corresponding with that of Vidal^^, as, 

• . . A^r^^A^rgA^rjA^r^-As . . . , (23) 

where K-i^n^ = A„, r% = B^„A-\, and 

Ab = A„_i. This is remarkable, because the iTEBD 
algorithm was originally conceived as evolving, simulta- 
neously, all bonds on an infinite-size translationally in- 
variant wavefunction, but it is formally equivalent to a 
lattice-growth scheme where two sites are added to the 
lattice at each iteration and the evolution is done per- 
formed via a single bond operator at the center of the 
lattice. This latter interpretation has advantages, in par- 
ticular the basis states in the iDMRG point of view are al- 
ways orthogonal, which provides a framework from which 
the convergence and stability properties are amenable to 
analysis. 

A. The orthogonality fideHty 

A measure of how close the wavefunction is to the fixed 
point is obtained by comparing the wavefunctions at suc- 
cessive iterations. At the fixed point, A^A~]_^ will be a 
diagonal matrix that represents the phases of the eigen- 
vectors chosen in the truncation procedure (if there are 
degenerate eigenvalues then the corresponding block can 
be an arbitrary unitary matrix) that can be incorporated 

into the neighboring matrices A^_^^ so that the wave- 
function transformation is the identity operation. Sim- 
ilarly, the diagonal matrix of phases of A~lj^A^ can be 
incorporated into -B^!!^Y- A measure for how similar the 
wavefunctions are can easily be obtained by calculating 
the fidelity between these wavefunctions at successive it- 
erations, which I call the orthogonality fidelity, 

i^ortho = (A„|A„_i) , (24) 

which is best obtained via the reduced density matrices. 




with the density matrices given by — A^^A^ and 
p„_i — AJj_j^A„_2- This reduces to the sum of the sin- 
gular values of A^A|^_^, so is quite easy and efficient to 
calculate. From the TEBD point of view, this measure 



shows how close the basis states are to being orthogonal. 
If there is any discarded weight in the truncation proce- 
dure of Fig. [T]step 5, this has a corresponding effect on 
-Fortho- That is, in iDMRG the wavefunction only reaches 
an exactly orthogonal fixed point if it is arranged that the 
truncation error is zero, or in iTEBD the requirement for 
exact orthogonality is that the time-step goes to zero. 

This points to an undesirable side-effect of a non- 
zero truncation of the wavefunction. While having a 
finite truncation is generally highly beneficial to con- 
vergence compared with, eg. single-site algorithms with 
no density-matrix mixing that are highly susceptible to 
convergence problems, near the fixed point of iDMRG 
a finite truncation error is a hindrance. To investigate 
the effects of truncation errors on the fixed point, I have 
compared iTEBD for various time-steps as well as two- 
site and single-site iDMRG algorithms with respect to 
how close the resulting state is to optimal, for a given 
number of states. Much of this section applies as well 
for finite-size calculations, although probably with less 
extreme differences. 

Compared with standard two-site (i)DMRG, the 
(i)TEBD algorithm with a careful reduction in time-step 
towards zero, produces a better wavefunction for a given 
number of states. Surprisingly, this difference is quite 
substantial, in all cases I looked at amounting to around 
30% improvement in the number of states needed for a 
given accuracy. In most applications of DMRG this is of 
no consequence, and it is not an efficient use of computer 
resources to try to converge to the best possible wave- 
function for a given number of states. Instead of perform- 
ing many many iterations to get 'optimal' convergence, 
often a better variational state is achieved with a less- well 
converged wavefunction by a small increase in the basis 
size. Nevertheless in some cases the number of states is 
all-important. In particular, to calculate an optimal state 
for a given m in a scaling analysis, an error bar of 30% is 
unacceptable. In addition, for applications such as time 
evolution, where the number of states kept grows sub- 
stantially (often exponentially) during the calculation, a 
30% reduction in the number of states required for the 
initial state may lead to a substantial improvement in 
the efficiency of the evolution. In a time-evolution calcu- 
lation, obtaining the initial state is usually much faster 
than the time evolution itself, so additional time spent 
on improving the initial state is worthwhile. 

Takasaki, Hikihara and Nishino^ have proposed a 
combination of two-site and single-site algorithms, us- 
ing the single-site algorithm in the last few sweeps of 
the calculation. However, this still doesn't produce the 
best possible result because the best states to keep do 
not correspond exactly with the largest m eigenstatcs of 
the density matrix. To be more precise, we can inter- 
pret the two-site algorithm as equivalent to a single-site 
algorithm with an environment that is perturbed by ad- 
ditional degrees of freedom corresponding to the second 
site. This perturbed environment means that the state 
isn't exactly translationally invariant due to the appear- 
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ance of a shallow bound-state at the central two sites, as 
noted by Dukelsky et al.^ But this perturbation to the 
wavefunction is such that the largest m density matrix 
states do not give the optimal states to keep at the point 
of crossing over to the single-site algorithm. With no 
good quantum numbers this is not such a serious prob- 
lem and I would expect that the algorithm converges to a 
nearly-optimal state regardless (although possibly quite 
slowly) . However practically all problems of interest uti- 
lize good quantum numbers, and in that case the choice 
of the number of states in each symmetry sector, which 
is fixed and does not change during the course of the 
single-site algorithm, will not be optimal. The solution 
is instead to evolve the wavefunction slowly, such that as 
the fixed point is approached the change in wavefunction 
at each iteration smoothly goes to zero. This does not 
happen in two-site DMRG because at the 'fixed point' the 
change in wavefunction per iteration is precisely the trun- 
cation error, which is usually small but non-negligible. A 
slow evolution of the wavefunction corresponds to imag- 
inary time-evolution with a small time-step, which could 
be done either with TEBD via the Trotter-Suzuki de- 
composition, or directly with DMRG by replacing the 
eigensolver by a multiplication with 1 — 5tH . A simi- 
lar effect could be obtained simply by performing just 
a single iteration of the eigensolver (Lanczos or similar) 
in the final sweeps of standard two-site DMRG. An al- 
ternative approach that should achieve the same effect 
is White's single-site density matrix mixing approacb^.^ 
and smoothly sending the density matrix perturbation 
to zero in the final stages of the calculation. 

B. Explicit orthogonalization 

The previous section has argued that, regardless of 
which algorithm is used to obtain the wavefunction, a 
finite truncation of the wavefunction leads a loss of or- 
thogonality, characterized by an orthogonality fidelity 
1 — -Fortho which will be of the same magnitude as the dis- 
carded weight of the truncation. This is not ideal, since 
even a small non-orthogonality affects the computation 
of expectation values, and the iterative algorithms for ex- 
pectation values described in Sec. IIVI below are simpler 
and perform better when the wavefunction is exactly or- 
thogonalized. Fortunately, this presents no problem as a 
procedure for explicitly orthogonalizing a translationally 
invariant MPS has been developed by Orus and Vidal^^. 

For concreteness let us take a unit cell of 2 sites, so the 
wavefunction at some iteration n is given by the matrices 

Ar^ KnB^ together with the center matrix at the previ- 
ous iteration, A„_i. Rotating the center matrix to the 
right gives the unit cell in the form of left-orthogonalized 
matrices, 

and the infinite wavefunction comprises this set of ma- 
trices repeated, ie, with P = A^A^lj^ the wavefunction 



is 

^n+l-^^n ^n+l-^^n ■ ■ ■ ■ ' ) 

If the wavefunction is exactly orthogonal, then P will be 
the identity matrix (up to a unitary matrix that could 
be incorporated into one of the A-matrices), and the or- 
thogonality fidehty will be 1 — Fortho = 0. But if P 7^ 1, 
the MPS is not orthogonal because no matter whether 
we attach P to the j4-matrix on the left or the right, 
the resulting A-matrix will not satisfy the orthogonal- 
ity constraint Eq. A simple-minded way of making 
the MPS orthogonal would be to simply rotate the P 
matrices to the right iteratively, until a fixed point of 
the updated ^-matrices is reached 4^ This is clearly not 
efficient, however. Instead, the procedure of Oriis and 
VidaP^ is to find a basis such that the dominant eigen- 
matrix of the transfer operator is the identity matrix*. 
The transfer operator in the left-orthogonal basis Tl is 
simply the identity operator of the unit cell, given by, 

Tl{E) ^ P^K%K'^EA^n'K\iP ■ (28) 

For an orthogonal state, Eq. ^ implies Tl(1) = 1 is 
a fixed point, as we would expect for a representation 
of the identity operator. For a non-orthogonal state we 
have instead Ti(l) ^ P^P. So, in order to transform 
to a basis in which P = 1, we obtain, with an efficient 
eigensolver such as the Arnoldi method'^'^, the dominant 
eigenvector Vl of T^, which we decompose into 

Vl ^ X^X . (29) 

Note that Vl is always Hermitian and non-negative. The 
factorization Eq. ([M)) could be done by a Cholesky fa- 
corization, however this is unstable if Vl is close to 
singular'^-; factorizing via a singular value or eigenvalue 
decomposition is more robust. We now perform a simi- 
larity transformation on the unit cell, 

A^l A^^^.P ^ XA^l Al\,PX-' , (30) 

and in this new basis, Pj,(l) = 1 as required. Note that to 
obtain the new A matrices, an additional step is needed. 
This could be done by solving the fixed point relation 
Eq. ^ but here we follow the outline of Ref. [H and 
solve the corresponding right-orthogonal transfer matrix 
for the right-orthogonalized unit cell, 

Tr{E) = QB'r.'Bl\,EBl%B^,,-^Q^ , (31) 

•Sl,S2 

where Q = K~^-yk^. Similarly to the left eigenvector, 
we now obtain the eigenvector Vr of Tr, and decompose 
into 

Vl = YY^ , (32) 
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from which we see that the complete set of transforma- 
tions is, 

A„_i ^ Xkn-lY (33) 

^ A^y 

n n 

A^ ^ XA^ 

This choice of tensors has the advantage over that used in 
Ref. m that the inverses of the X, Y matrices are not ac- 
tually required, only the inverse of XA„_iF which deter- 
mines the singular values in the orthogonal basis. There 
is one other difference here to the procedure described by 
Orus and Vidal, concerning the treatment of unit cells of 
> 1 sites. Orus and Vidal suggested handling these 
unit cells with a coarse-graining procedure that combines 
a unit cell of N cites into a single site with an increased 
local dimension. This is inefficient, as it corresponds to 
evaluating the summations in the eigenvalue equation Eq. 
(|28p in a non-optimal way, using 0{d^) matrix multi- 
plies rather than 0{Nd) multiplies if the summations 
are performed sequentially. Aside from this difference, 
and the slightly different choice of tensors, the proce- 
dure described here is equivalent to that of Orus and 
Vidal. Note that this procedure orthogonalizes the unit 
cell, it does not however guarantee that the individual 
sites comprising a multi-site unit cell separately satisfy 
Eq. (2]). The step of orthogonalizing each site should be 
done for the sake of numerical robustness, and is readily 
achieved through the usual orthogonalization procedure, 
described in Sec. HIl 

The orthogonalization procedure is quite efficient, and 
for robustness I always perform this orthogonalization 
step before calculating the final expectation values. 

IV. EXPECTATION VALUES 

In this section, we consider expectation values of 
MPO's, firstly for the case of an operator that acts only 
on a finite range of the lattice, then for translationally 
invariant triangular MPO's. 

A. Finite-range operators 

Any finite-range operator of width N can be expressed 
as a position-dependent MPO with matrices Wi . . . Wn 
with an open-boundary condition representation. That 
is, the left-hand MPO matrix is a 1 x M row vector, and 
the right-hand MPO is a M' x 1 column vector. Starting 
from the right-hand side, the initial operator represents 
the identity operator on an infinite lattice, which will 
be trivial if the state is orthonormalized. Thus define 
Fq = 1 to be this identity operator, then starting from 
the right hand side, 

Fn+i = E W(^%i)o.'.B'' KB'^ . (34) 



The final iteration results in a single matrix that rep- 
resents the operator matrix elements. The expectation 
value is then given by TrF^p as usual, where p is the 
density matrix for the semi-infinite strip, being the fixed 
point of p = Y^sB'^pB'. 

B. Triangular operators 

To evaluate an expectation value of a lower-triangular 
MPO, we need to determine the matrix elements as a 
function of the wavefunction A-matrix elements. For an 
MPO of dimension M x M, we have M different matri- 
ces i?", a — 1,2,..., M, that represent the components 
of the MPO on the remainder of the lattice. For con- 
creteness, we firstly assume that the operator is of the 
form of a sum of short-range terms. This means that 
the diagonal of the MPO is zero, except for the top-left 
and bottom-right entries which are equal to the identity 
operator. This gives the fixed-point recursion 

E^' = E wil^A^'^E^A' - EoSc.,,1 , (35) 

where Eq is the expectation value per site. This is 
straightforward to evaluate term by term. Starting from 
the last entry at a = M, since the bottom-right entry is 
the local identity operator we have, 

E'^' = E ^"''E'^'A' , (36) 

s 

which is just the fixed point of the orthogonalization con- 
straint, so for left-normalized matrices E^^ = 1 identi- 
cally. The remaining matrices E" for 1 < a < M are 
simple to construct in descending order because these 
matrices are functions only of the previous matrices i?" 
for a' > a. The exception is the final matrix E^, which 
again has an identity operator on the diagonal, giving the 
linear equations 

E^-^A''^E^A' + Eo^ J2 wfjA''lE°'A' , {37) 

s a>l,s' ,s 

where the expectation value Eq is given by the trace of 
the right hand side. If only the expectation value is re- 
quired, these linear equations do not need to be actually 
solved. The full solution is required only if the matrix 
E^ needed, for example to restart an iDMRG calculation 
from a given translationally invariant wavefunction. 

The procedure when there are additional matrix ele- 
ments on the diagonal of the MPO follows similarly. 

C. Transfer matrix and symmetry breaking 

The spectrum of the transfer matrix Eq. (|28p is very 
useful as this determines the possible correlation lengths 
in the system.^ The transfer matrix has always at least 
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one eigenvalue equal to 1, as follows from the normal- 
ization constraint Eq. The next leading eigenvalue 
determines the longest correlation length in the system. 
Normally this eigenvalue will be strictly less than 1. Even 
in the presence of long-range order, a non-zero order pa- 
rameter will lead to a spontaneously broken symmetry 
at the infinite-size fixed point, as long as the symme- 
tries that are imposed on the MPS representation allow 
this. In an infinite size system, a broken symmetry im- 
plies multiple groundstates which are exactly degenerate. 
Such spontaneous symmetry breaking typically does not 
occur in finite-sized systems because there are usually 
matrix elements connecting the groundstates that lift the 
degeneracy, and such matrix elements are exponentially 
small in the lattice size ~ exp[— L]. Note that trans- 
lation symmetry is an exception to this, since a finite 
system inherently breaks translation symmetry from the 
beginning. 

Spontaneous symmetry breaking of discrete symme- 
tries is realized in an infinite MPS with matrix elements 
connecting degenerate groundstates going to exactly zero 
at the fixed point. Hence we expect in a numerical sim- 
ulation that exactly one of the degenerate states will be 
chosen non-deterministically. (If any of the groundstates 
has a lower entanglement then this state would tend to 
be preferentially selected by the variational principle, 
as this state would have a lower energy for a particu- 
lar basis size.) However, if good quantum numbers are 
used then the obtained groundstate will always have this 
symmetry even if it would otherwise be spontaneously 
broken. For example. In the low-field ordered phase of 
the ITF model Eq. ([H]) . there are two degenerate ther- 
modynamic groundstates, namely the two fully polar- 
ized states I ii . . .} and | f j . . .). In the absence of any 
mechanism enforcing Z2 symmetry, iTEBD or iDMRG 
will spontaneously select one of these groundstates, se- 
lected eg. by numerical noise from finite-precision arith- 
metic that leads to one groundstate having slightly higher 
weight than the other, which then quickly dominates.'*^ 
However if a Z2 quantum number is used, then the ob- 
tained groundstate will be an exact eigenstate of Z2, 
namely | i| . . .} ± | H . . .), which is a GHZ stated with 
long-range order leading to an additional eigenvalue of 
unity in the transfer matrix spectrum. Thus, the transfer 
matrix is very useful for diagnosing the presence of long- 
range order, without a priori knowledge of the specific or- 
der parameter. Similar considerations also apply for the 
transfer matrix extended to string order parameterS)^^ 
where the spontaneous symmetry breaking is instead as- 
sociated with edge excitations that do not affect the bulk 
wavefunction^S. 

Figure[3]shows the transfer matrix spectrum of the ITF 
model in the vicinity of the critical point A = 1, for some 
different basis sizes. As can be seen, the ITF model, 
which is remarkable for having a very large correlation 
length even with a very small basis size, has relatively 
few eigenvalues that contribute to long range correla- 
tions. For example, considering only the eigenvalues cor- 




FIG. 4: The spectrum of the transfer matrix for the ITF 
model near the critical point. The next-leading eigenvalue de- 
termines the correlation length, here ^ — 578.7, 1236.6, 2539.2 
for m = 20, 30, 40 respectively. 



responding to correlation lengths of greater than 10 sites 
gives only around 30 relevant eigenvalues, even for the 
case of TO = 40 where the correlation length is a few thou- 
sand sites. Thus the long-distance behavior of correlation 
functions can be determined by looking at just a few 
eigenvalues. In the case of a critical system, the power- 
law behavior of a correlation function is obtained approx- 
imately, as the sum of a set of exponentials^. Hence the 
eigenvalues of the transfer matrix give a good guide as 
to the range of validity where an approximate power- 
law can be expected, ie. to a distance the maximum 
correlation length in the system. A complication of us- 
ing the transfer matrix directly to expand a correlation 
function as a set of exponentials is that the transfer ma- 
trix of Eq. (|28[) is not a normal matrix and hence the 
eigenvectors are not orthogonal. But this presents no es- 
sential difficulty, as long as care is taken to project out 
the large number of components with short correlation 
lengths prior to determining the amplitudes (expansion 
coefficients) of each mode, eg. by calculating the operator 
i?-matrix for a correlation at some distance larger than 
the smallest correlation length eigenvalue obtained from 
the transfer operator. 



D. Fidelity 

The overlap of two MPS is a useful operation that 
is used to obtain the groundstate fidelity F{gi,g2) = 
I (gi 152) I between groundstates of different interaction pa- 
rameters gi and §2- The fidelity has attracted much 
recent attention as a tool for probing quantum phase 
transitions*^^ For an MPS, this quantity is always of the 
form d^, for an L-site lattice with < d < 1. Thus, 
the fidelity is measured per site in much the same way 
as the expectation value of an exponential MPO (indeed, 
the fidelity is equivalent to the expectation value of the 
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identity operator which can be seen as the simplest pos- 
sible exponential MPO). The fidelity susceptibility'^^ is a 
measure of the curvature of the fidelity at a particular 
point, and hence shows how quickly the wavefunction is 
changing with the interaction parameter. The fidelity 
susceptibility per site is given by, 



XF{g) = 



lim 



-21nf(ff,g + ^g) 
LSg^ 



(38) 



For finite MPS, the fidelity has been proposed as a mea- 
sure since at least the work of Andersson, Boman and 
Ostlundf2S For an infinite MPS, the fidelity (and hence 
also the fidelity susceptibility) is straightforward to cal- 
culate via an eigenvalue equation. Similar to the transfer 
operator (which could be seen as calculating the norm 
(^fl^)), the fidehty between two states | ^1) and | *2) is 
given by the largest eigenvalue of an overlap operator, 
here given for a two-site unit cell example (compare Eq. 
|), assuming the states are orthogonalized) , 

T,,{E) = J2 ^r^(f)Ar^(l)i?Ar(2)A^^(2) , (39) 



where A*"(l) and ^^^"(2) are the matrices representing 
I'^l) and |^'2) respectively. The fidelity per unit cell is 
then given by the eigenvalue A of Ti2{E) — XE, which 
can be any complex number of magnitude < 1. 

In a finite MPS, the fidelity is sometimes complicated 
by boundary effects that are irrelevant to the bulk prop- 
erties of the wavefunction. For a concrete example, 
the bilinear-biquadratic S = 1 chain*^"^ has a variety of 
groundstate phases including a gapped resonating va- 
lence bond (RVB) phase (including the AKLT point"*^ 
that is exactly solved with a non-Abelian MPS with di- 
mension m = 1), a gapped dimer phase, and a critical 
trimerized phase. The RVB phase is characterized by 
effective spin- 1/2 edge excitations that are best treated 
by applying real spin-1/2 states at the boundaries, which 
means that the bulk wavefunction in the infinite size limit 
has quantum numbers that are all half-integer. This 
presents no problems for MPS algorithms in the infi- 
nite limit, however in other phases such as the gapped 
dimer phase it is best not to supply spin-1/2 states at 
the boundaries, as these frustrate the formation of dimer 
pairs and also introduces spurious effects into the fidelity. 
On the other hand, to calculate the fidelity for a finite 
size system the Hilbert space of the wavefunctions must 
coincide exactly, which causes severe difficulties in calcu- 
lating the fidelity across the RVB-dimer phase transition. 
These problems do not arise in the infinite MPS case, as 
the boundary effects are not present and having different 
spin representations occurring in the bond basis of the 
wavefunctions is no problem for the fidelity calculation. 



V. CONCLUSIONS 

In this paper I have examined the infinite-size DMRG 
algorithm from the perspective of the matrix prod- 



uct approach, which gives directly a translationally in- 
variant fixed point previously studied by Rommer and 
Ostlund^i^. The core development in this algorithm is 
the wavefunction transformation, which improves upon 
that used in PWFRG^'^, and is useful even for small fi- 
nite systems. Sec. llVl has shown how to calculate expec- 
tation values in the thermodynamic limit of operators 
directly from the matrix product form of an observable, 
which gives access to higher moments such as the energy 
variance. The spectrum of the transfer matrix gives de- 
tailed information on the correlation length and ordered 
phases, and expanding a correlation function as a sum of 
exponentials corresponding to the modes of the transfer 
matrix gives an alternative approach to fitting a criti- 
cal exponent. The groundstate fidelity between different 
groundstates is readily calculated, and this gives a precise 
way to determine the fidelity that is free of complications 
due to boundary effects. 

Infinite-size DMRG represents an efficient method for 
finding translationally invariant states that could replace 
many of the current uses for finite-size DMRG. How- 
ever one situation in which the current algorithms are 
not a sufficient replacement is the calculation of excited 
states, for example energy gaps, low-lying excited states 
and spectral functions. Despite some suggestions to the 
contrary^^, I believe these calculations should be readily 
achievable. Indeed, the basic approach of representing a 
particle excitation as a superposition of local perturba- 
tions with an MPS has already been discussed by Rom- 
mer and Ostlund^, leading to an ansatz for an excitation, 
for example the application of a particle creation cj. at 
some momentum k to an infinite MPS. Building upon 
this approach I see no reason why there should not exist 
efficient infinite-size counterparts for the well-established 
correction vector-^ and dynamical DMRG'^'^ algorithms 
for frequency-space dynamical correlations of finite-size 
systems. Note that calculating spectral functions via 
real-time evolution of a local particle excitationi^ is al- 
ready achievable, by using a section of a translationally 
invariant state as the input for a finite-size time evolu- 
tion. 

The iDMRG algorithm presents an interesting di- 
chotomy of solving simultaneously for both the wavefunc- 
tion A-matrix elements and the Hamiltonian iJ-matrix 
elements. The present scheme uses an efficient Lanczos 
eigensolver to obtain the wavefunction, while the effect 
of growing the lattice one site at a time amounts to using 
the less efficient power method to obtain the Hamiltonian 
matrix elements, in the sense that exactly one iteration of 
the _B-matrix multiply is performed per step. Therefore, 
this algorithm will converge in a minimum number of it- 
erations that is at least as large as the correlation length. 
On the other hand. Sec. IIVBI shows how to obtain the 
Hamiltonian matrix elements of a translationally invari- 
ant state much more efficiently using a linear solver. This 
suggests an accelerated algorithm that efficiently solves 
simultaneously the fixed point wavefunction and the as- 
sociated fixed point Hamiltonian Eq. (|35p , which could in 
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principle converge in a number of iterations much smaller 
than the correlation length at the cost of the algorithm 
being non-linear. 

Algorithms for finding the groundstate of a one- 
dimensional Hamiltonian form a component of some al- 
gorithms for higher dimensional systems, most notably 
the iPEPS algorithm-4 for 2D tensor product states^. 
The higher efhciency of direct eigensolvers should serve 
a similar function as the ID case, although the analogue 
of a triangular MPO for the Hamiltonian construction is 



clearly more complicated in 2Di^ 
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